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ABSTRACT 

In this pape r we extend our numerica l meth od for simulating terrestrial planet for- 
mation from Leinhardt & Richardson (2005) to include dynamical friction from the 
unresolved debris component. In the previous work we implemented a rubble pile plan- 
etesimal collision model into direct iV-body simulations of terrestrial planet formation. 
The new collision model treated both accretion and erosion of planetesimals but did 
not include dynamical friction from debris particles smaller than the resolution limit 
for the simulation. By extending our numerical model to include dynamical friction 
from the unresolved debris, we can simulate the dynamical effect of debris produced 
during collisions and can also investigate the effect of initial debris mass on terres- 
trial planet formation. We find that significant initial debris mass, 10% or more of 
the total disk mass, changes the mode of planetesimal growth. Specifically, planetesi- 
mals in this situation do not go through a runaway growth phase. Instead they grow 
concurrently, similar to oligarchic growth. The dynamical friction from the unresolved 
debris damps the eccentricities of the planetesimals, reducing the mean impact speeds 
and causing all collisions to result in merging with no mass loss. As a result there 
is no debris production. The mass in debris slowly decreases with time. In addition 
to including the dynamical friction from the unresolved debris, we have implemented 
particle tracking as a proxy for monitoring compositional mixing. Although there is 
much less mixing due to collisions and gravitational scattering when dynamical fric- 
tion of the background debris is included, there is significant inward migration of the 
largest protoplanets in the most extreme initial conditions (for which the initial mass 
in unresolved debris is at least equal to the mass in resolved planetesimals). 



1 INTRODUCTION 



Extrasolar planets are numerous and diverse, with recently 
detected examples ranging from "hot Jupiters", with peri- 
ods of days, to super-Earths, with masses 5M^ (Anderson 
et al.|2008| [Mayor et al.|2009[ ). It is evident that planet for- 
mation is common and ubiquitous. However, results from 
numerical simulations lag behind recent observational dis- 
coveries. Simulations of the formation of our own solar sys- 
tem often result in planetary eccentricities that are too high 



and/or ejection of one of the terrestrial planets (Raymond 
et al. 20081. Damping of eccentricities requires either gas 



is to begin with a simple model and systematically make 
the model more and more realistic. In this paper, as with 
our previous paper (Leinhardt & Richardson|2005|, we have 
chosen to focus on the terrestrial region during the middle 
phase of planet formation, assuming no gas but allowing for 
the production of background debris. Thus in this work we 
ignore the effects of gas. In this context we use an A^-body 
code to model the coUisional and dynamical evolution. 



1.1 Previous Work 



( Tanaka fc W ard 2004) or large amounts of small debris iKokubo & Idal ([2002| completed a series of direct A-body 



( Goldreich et al.(2004| ). It has yet to be shown quantitatively 
whether either situation arose in our own solar system. 

Numerical simulations permit modeling the complex in- 
terplay of physical processes during planet formation that 
are otherwise difficult to assess. Observations are snapshots 
that provide only limited information on how a protoplan- 
etary disk evolves into a solar system. In addition, not all 
phases of planet formation are observable. Our approach to 
understanding the formation and evolution of solar systems 



simulations of the middle phase of terrestrial planet forma- 
tion (for reviews see |Lissauer|1993]|Chambers|2004| ). In these 
simulations it was assumed that every collision resulted in 
perfect merging (no mass loss), thus, every collision resulted 
in growth of the planetesimals. There was no erosion nor the 
possibility of producing debris that could damp the larger 
protoplanets found at the end of the simulations. In our 
previous paper we developed a more realistic planetesimal 
collision model that allowed both accretion and erosion. We 
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then ran evolution models with similar starting conditions 
to Kokubo & Ida (20021. We found virtually identical re- 



dius with bulk density 2 g cm ) . As a result, these planetes- 
imals are in the gravity regime where their material strength 



suits. The number, mass, and spatial distribution of the 
protoplanets were similar and showed both runaway and oli- 
garchic growth. However, this earlier work did not include 
a model for the full dynamical feedback of the coUisional 
debris (mass elements below the simulation resolution) on 
the planetesimals. Though planetesimals could accrete coUi- 
sional debris, there was no dynamical friction from the debris 
on the planetesimals. In addition, all simulations in |Lein-| 
hardt & Richardson ( 2005 1 used a low mass of initial debris: 



the effect of starting with a larger amount of mass as debris 
was not investigated. 

The middle phase of planet formation is not observable 
(planetesimals are too big to be observed by infrared obser- 
vations and too small to be observed at visible wavelengths). 
In addition, the physical mechanisms that govern the ear- 
lier phases, namely planetesimal formation, are disputed (for 
example, the long-standing debate between gravitational in- 
stability and turbulence models; [Johansen et al.|2007l|CuzzIl 
et al. 20081. As a result, the initial conditions for numeri- 



cal simulations of this phase are unknown. In general, it is 
usually assumed that planetesimals that have just decou- 
pled from the gas are on similar, almost circular orbits with 
small inclinations. However, the amount of material in plan- 
etesimals and the amount in smaller debris is completely 
unknown, and as with many of the protoplanetary disk pa- 
rameters, the amount of debris may vary from star to star 



depending on metallicity and/or spectral type. Leinhardt 
& Richardson ( 2005| l chose an initially low mass in debris 



because they did not have full feedback in their numerical 
model. 

In this paper we investigate the effect of the initial de- 
bris distribution on the middle phase of planet formation. 
We use a similar numerical method as Leinhardt & Richard- 



( 2005 1 but one that is upgraded to include dynamical 



friction from unresolved material and the ability to track 
compositional mixing in both the planetesimals and the de- 
bris through the course of the simulation. 



2 METHOD 

For our simulations we use the parallel A^'-body gravity code 
pkdgrav. The code uses a second-order leap-frog integrator 
with a hierarchical tree ( Richardson et al.|2000 Stadel|2001 1 
to provide computation time scaling as 0{N log N), where 
N is the number of resolved planetesimals. The version used 
here has been modified to include realistic collisions between 



planetesimals (see |Leinhardt fc Ric hards on 2005) for details) 
and accounts for dynamical friction from the collisional de- 
bris. 



2.1 Planetesimal Model 

In this section we summarize the planetesimal structure and 
collision model used in the numerical simulations. The plan- 
etesimal model used in this paper is the same as that used in 



Leinhardt & Richardson (20051. Please refer to that paper 



for more details. 

Due to numerical limitations, our direct numerical sim- 
ulations must start with large planetesimals (-^ 60 km in ra- 



is negligible compared to their gravitational strength (Lein- 
hardt et al.[2008 , Leinhardt fc Stewart|2009| [Stewart fc Leii> 



hardt 2009 and references therein). Thus, we have chosen 
to model the planetesimals in our numerical simulations as 
gravitational aggregates (or, more precisely, idealized rubble 
piles — cf. Richardson et al.|2002 | during collisions. 

In the numerical simulations presented here, plan- 
etesimals grow via accretion of debris and planetesimal- 
planetesimal collisions. The collisions between planetesimals 
are compared to a look-up table of collision outcomes, and 
either the tabulated result is used (if suitable), or a direct 
simulation of the interaction is performed. Collision param- 
eters of impact speed, impact parameter, and mass ratio, 
as well as a user-specified coefficient of restitution are used 
to interpolate or extrapolate the collision outcome from the 
collision outcome database. The database contains the mass 
of the largest post-collision remnant from several hundred 
rubble-pile planetesimal collision simulations over a wide 
range of parameter space ( Leinhardt et al.]|2000 Leinhardt 
|fc Richardson,2002_, 2005| . 

If the predicted collision outcome from the database is 
one large remnant and a small amount of debris, the colliding 
particles are replaced by the large remnant and the debris is 
followed semi-analytically (see below) . If, on the other hand, 
the collision results in two massive bodies (the largest post- 
collision remnant being less than 80% and second largest 
remnant being greater than 20% of the total mass of the 
system), the planetesimals — which have been modeled as 
single particles up to this point — are replaced by rubble piles 
of 100 or so particles each, and the collision is integrated 
explicitly in the planetesimal disk. 

In this latter case of a resolved collision, for the first 
ten dynamical times (Tdyn ~ where G is the gravi- 

tational constant and p is the bulk density of the planetesi- 
mals), the particles are only allowed to bounce inelastically 
(there is no merging, in order to allow the system to relax in 
a realistic way, nor fracturing, because the re-impact speeds 
between particles are small). After ten dynamical times, the 
rubble particles merge with one another when they collide. 
After twenty dynamical times, any particle from the colli- 
sion that is smaller than the resolution limit (the size of a 
planetesimal at the start of the simulation) is demoted to 
unresolved debris. The time allotted for these bouncing and 
merging phases is based on studying the collision experi- 
ments used to generate the outcome database. It was found 
that in most cases intervals of ten dynamical times were 
an adequate compromise between detailed collision outcome 
modeling and computational expediency. 

In our model, the planetesimal disk is divided into a 
configurable number of annuli for the purpose of tracking 
debris. Unresolved debris from a collision is added to the 
mass density of the annulus at the location where the debris 
was generated. The debris is assumed to have circular Ke- 
plerian orbits and a fixed scale height of 1 x 10"'' AU. The 
value of the scale height was chosen such that all coplanar 
planetesimals are fully embedded in the debris disk during 
the entire evolution of the system. The planetesimals sweep 



up debris as they pass through the annuli (see §2.3 of Lein- 
[hardt fc Richardson|2005l) and are damped by the dynamical 
friction of the debris {[ 2.2 below). If a collision occurs at the 
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edge of the debris simulation domain (either interior to the 
inner edge of the inner annulus or exterior to the outer edge 
of the outer annulus) any unresolved debris outside the de- 
bris simulation domain after twenty dynamical times is put 
into a "trash bin" for tracking purposes. The mass density 
and composition of the trash bin are tracked through the 
entire simulation like any other annulus bin but the con- 
tents do not interact with the planetesimals. The resolved 
planetesimals are free to move inside or outside of the de- 
bris simulation domain. If they move outside, the resolved 
planetesimals will not interact with any debris during that 
interval. 



In the simulations presented in this paper we have sped 
up the planetesimal growth by artificially inflating the radii 



of the planetesimals by an expansion parameter / | Kokubo 
ifc Ida_2002; Leinhardt fc Richardson|20b5| ). As a result, the 
dynamical friction (Eq. |4|) must be modified as well in or- 
der to keep pace with the accelerated collisional evolution. 
The collision cross-section scales roughly as the surface area 
for most of the evolution modeled in these simulations until 
gravitational focusing becomes important. Thus, the dynam- 
ical friction impulse applied becomes 



AVM = ^ q / AiVM 



(5) 



2.2 Dynamical Friction 

Unresolved debris circularizes the planetesimals through the 
action of periodic impulses via equation 7-14 of |Binney fc| 
Tremainel (119871), 



dt 



-vm(1) 



where the vm is the difference between the orbital velocity 
of the resolved planetesimal and that of the unresolved de- 
bris at the planetesimal's instantaneous position, m is the 
mass of a single unresolved debris particle, M is the mass 
of the planetesimal, and the integral gives the amount of 
unresolved material with speed less than the speed of the 
planetesimal, for an assumed debris particle speed distribu- 
tion f{Vm)- In the Coulomb logarithm for a Keplerian disk, 
A = bmax{vl,i + 2Vdispvlirc)/iGM), the maximum impact 
parameter between a planetesimal and the unresolved debris 
is bmax = Rh + ay/Vdisp + siu^ i, whcrc Rh = isM^y^^ is 
the Hill radius of the planetesimal, a is the semi-major axis 
of the planetesimal, i is the inclination of the planetesimal, 
Vdisp is the velocity dispersion of the unresolved debris, Vdrc 
is the instantaneous circular speed of the planetesimal, and 



G is the gravitational constant (Stewart & Ida 2000 Ford 
fc Chiang|20"07l ) 



In this paper we assume the unresolved debris is "cold" 
(having zero eccentricity and inclination). Thus, the relative 
speed of the planetesimal to the mean Keplerian speed of 
the debris disk is generally much larger than the velocity 
dispersion of the unresolved debris. Let us assume that the 
unresolved debris has a Maxwellian distribution (Eq. 7-16 
Binney fc Tremaine|1987| , 

no 

— exr> I — ^ . 



exp(- 



1 2 /■,,2 X 



(2) 



where no is the number density of unresolved debris parti- 
cles. We may then integrate Eq.[l](see equation 7-17 Binney 
[&Trcmainc 1987), 



-2Trln{l + A:^)G^pM 



dt 



erf(X) 



|vm, (3) 



where X = vm /{V2Vdisp), P ~ nom is the mass density 
of the unresolved debris particles, and m <C M, therefore, 
M+m ~ M. We now apply the assumption that the Vdiap ^ 
Vm (see !|4]for discussion of the relaxation of this condition), 
so X goes to infinity and Eq. [S] becomes, 

dvM _ -27rln(l-hA^)GVA/ 



dt 



-Vm- 



(4) 



where At is the time between impulse applications. In all 
simulations presented here, f — 6. At varies between 1 and 
10 timesteps depending on the initial mass of the debris and 
the magnitude of the resulting impulse (see j |2.5| fc Table [l] 
for more details about time steps). 



2.3 Composition Tracking 

In addition to the implementation of full dynamical friction, 
we have added the capability to track mixing both in plan- 
etesimals and the unresolved debris. Each planetesimal car- 
ries a mass-weighted fractional composition histogram that 
is binned by the semi-major axis. The number of histogram 
bins is user-specified — in all simulations presented here we 
used 15 bins. Each annulus of unresolved debris carries an 
analogous composition histogram. At time zero, all bins in 
each histogram have zero value except for the bin corre- 
sponding to the initial particle/debris location, which has 
a value of one. As a planetesimal grows, the populations of 
its histogram bins vary according to the evolving compo- 
sition. For debris, the histogram bins evolve as new debris 
enters into the annulus, or existing debris is swept up by 
planetesimals. 

Figure [l] shows a cartoon of a collision between two 
planetesimals (A fc B) in a swarm of background debris 
(small unlabeled circles). The composition histograms for 
the planetesimals and debris are shown below the cartoon 
of the collision. In the example, there are five semi-major 
axis bins. Planetesimal A is originally made up completely 
of material from the middle of the protoplanetary disk (blue 
rectangle in lower left of Fig.[T]). Planetesimal B is originally 
made up of material from the inner annulus (red rectangle), 
and the debris is originally made up entirely of material in 
between A and B (green rectangle). 

After the collision, there is one planetesimal, C, and 
some additional unresolved debris. Planetesimal C is a mass- 
weighted compositional mix of planetesimals A fc B. The 
composition histogram for C (lower right) has some mate- 
rial from planetesimal A (blue rectangle in the C row) and 
some from planetesimal B (red rectangle in the C row). The 
debris is also mixed as a result of the collision, now con- 
taining material from both original planetesimals (blue and 
red rectangles) in addition to the original debris material 
composition. 

A planetesimal that passes through the now heteroge- 
neous mixed debris will sweep some of it up. The composi- 
tion of the debris will modify the composition histogram of 
the growing planetesimal, just as if it had been modified in 
a collision with another planetesimal. 
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Figure 1. Cartoon of a collision between two planetesimals (top left), A (blue) and B (red), in a swarm of unresolved debris (smaller 
green circles). The result (top right) is planetesimal C (purple) and heterogeneous debris (small olive-color circles). Example composition 
histograms pre- and post-collision are shown on the bottom left and right, respectively. 



2.4 Planetesimal Disk Model 

In this paper we present moderate-resolution (initial num- 
ber of resolved planetesimals — 10*) simulations of var- 
ious initial debris disk masses to investigate the effect of 
environment on protoplanet formation (see ^3.11. For all 
simulations we used the standard model for the resolved 
planetesimal disk from Leinhardt & Richardson (20051. If 



all of the mass is in resolved planetesimals (Mr) and no 
mass is in unresolved debris (Md), the standard model re- 
duces to a "minimum- mass solar nebula". For all simula- 
tions, the initial surface density of planetesimals at 1 AU, 
El ~ 10.0 g cm~'^, and the surface density distribution 
Ep = Ei(a/1 AU)^^''^ The total mass of unresolved de- 
bris varied from — 25 , with surface density at 1 AU 
(E'l) between and 100 g cm~'^, and a surface density dis- 
tribution set to T^debris = E'l (o/ 1 AU) "^''^ . Each simulation 
began with a 1 AU-wide band of particles and unresolved 
debris centered at 1 AU. The simulations were run for at 
least 1 X 10^ yr — long enough to initiate runaway growth 
and see the formation of multiple protoplanets in previous 
work ( [Leinhardt fc Richardsonl|2005| [Kokubo fc Ida|2002 L 



2.5 Time Step 

In the numerical simulations presented here there are 
two time scales: the orbital time scale of the planetesi- 
mals (planetesimal-Sun interactions) and the time scale of 
planetesimal-planetesimal interactions. The orbital dynami- 



cal time is one year at 1 AU. The planetesimal-planetesimal 
interaction time scale is significantly less (~ 1 hr for p — 2 g 
cm""^) as it depends inversely on the bulk density of the col- 
liding planetesimals ( ^2.1[ |. The minimum step size is a func- 
tion of the choice of integrator and the dynamical time. For 
leapfrog integrations of low-eccentricity orbits, ~ 33 steps 
per orbit are required to accurately resolve the orbit (IQuinnl 



et al. 19971. In Leinhardt & Richardson (20051, we were 



more conservative and used ~ 100 steps per dynamical time 
to resolve orbital or collisional interactions between plan- 
etesimals. The large difference in the time step required to 
model collisions and orbits has led us to implement a bi- 
modal particle time stepping scheme. The large step is used 
to integrate the orbits of the planetesimals while the small 
steps are used during the interval between when a collision 



is detected (for details see §2.5 of Leinhardt & Richardson 



2005) and finally resolved. In this paper we have found that 
the addition of dynamical friction reduces the minimum time 
step required to accurately integrate the planetesimal orbits. 
In other words, 100 steps per orbit at 1 AU is not sufficient 
resolution for all initial conditions considered here (discussed 
below). This is because the impulse (Eq.[5| from the dynam- 
ical friction of the debris can be large in comparison to the 
speed of a given planetesimal, resulting in a large change 
in direction and/or magnitude of the planesimal's velocity 
vector. The step size required is strongly dependent on the 
magnitude of the dynamical friction impulse and thus the 
debris mass. 

For each set of initial conditions we completed a time 
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Table 1. Summary of simulation parameters. 



Time [yr] 

Figure 2. Number of particles versus time using four different 
time steps. Black: At = 2.5 X IQ—* yr, blue: At = 1 X 10""' yr, 
red: At = 0.5 X 10"'' yr, green: At = 0.25 X IQ— * yr. The curve 
for the simulation with the longest time step gradually separates 
from those using the smaller time steps, indicating that 2.5 X 10"'* 
yr is too large a step to accurately model the dynamical evolution 
of the protoplanetary system. 



step test that consisted of comparing the number of parti- 
cles as a function of time for simulations of the same initial 
condition but different major (i.e. orbital) time steps. The 
time step necessary for a particular initial condition was de- 
termined when the number of particles versus time agreed 
with a simulation of smaller time step. Figure |2] shows the 
number of particles versus time for a simulation with the ini- 
tial mass of the unresolved debris equal to 10% of the mass 
of planetesimals. The black, blue, red, and green points cor- 
respond to major steps of 2.5, 1, 0.5, and 0.25 x 10~* yr, 
respectively. The evolution of the number of particles in the 
largest step size simulation is slower than for the smaller 
step sizes. However, all three smaller step sizes have a simi- 
lar slope, thus for this simulation a time step of 1 x 10~* yr 
was used for the major step (see Table [l] for time steps used 
in each simulation). 

In addition to the time step test, the angular momen- 
tum and energy conservation was checked for a perfect merg- 
ing simulation that does not include the dissipative effects 
of dynamical friction from unresolved debris. The fractional 
angular momentum after 1 x 10^ yr was conserved to one 
part in 10^ and the energy to one part in 10*. When unre- 
solved debris is included, the angular momentum of the re- 
solved planetesimals increases as the planetesimals accrete 
unresolved debris, as expected. 



3 RESULTS 

Table[T] summarizes the chosen parameters of individual sim- 
ulations. In the table, the first column gives the simulation 
name while the second column indicates whether dynami- 



Sim. 


Dyn. Fric. 




At (yr) 


St (yr) 




Ttot (yr) 


la 


N 


0.00 


1 X 10-3 


1.25 X 10" 


4 


2.1 X 10^ 


lb 


Y 


0.00 


1 X 10-"* 


5.00 X 10- 


5 


1.0 X 10^ 


2a 


N 


1.00 


1 X 10-3 


1.25 X 10" 


4 


4.0 X 10^ 


2b 


Y 


1.00 


1 X 10-4 


5.00 X 10- 


5 


1.0 X 10^ 


3a 


N 


10.0 


1 X 10-"* 


5.00 X 10" 


5 


1.0 X 10^ 


3b 


Y 


10.0 


1 X 10-"* 


5.00 X 10- 


5 


1.0 X 10^ 


4a 


N 


100. 


1 X 10-* 


5.00 X 10" 


5 


1.0 X 10^ 


4b 


Y 


100. 


1 X 10-1 


5.00 X 10- 


5 


1.0 X 10^ 



cal friction from the unresolved debris was included or not 
([Y]es or [N]o). The rest of the column labels are as follows: 
is the surface mass density of the unresolved debris at 1 
AU, At is the major time step, St is the minor time step, 
and Ttot is the total simulation time. The control simula- 
tions (sims. l-4a) did not include dynamical friction from 
the unresolved debris and were run using the same numeri- 
cal method as |Leinhardt fc Richardson] ( |2005"| ). Comparisons 
between the control simulations and those with dynamical 
friction (l-4b) are presented in j 3.1|3.5 In [3.1 



we present 

a summary of the changes in planetesimal evolution that 
occur when dynamical friction of the unresolved debris is 
included in the model. In 93.21 we discuss the differences in 



the composition of the protoplanets after 10^ yr. The evolu- 



tion of the unresolved debris is presented in |3.3[ In §3.4| we 
show that significant initial debris mass changes the growth 
mode of the planetesimals when the dynamical friction of 



the unresolved debris is included. In the last section {i 3.5 1 
we discuss the differences in the collision outcomes between 
the control simulations and those that include dynamical 
friction of the unresolved debris. 



3.1 Effect of Dynamical Friction 

Dynamical friction from unresolved debris changes the evo- 
lutionary growth of protoplanets by damping eccentricities 
and, in some cases, causing significant inward migration. 
Figures |3] and |4] show the eccentricity and mass of the plan- 
etesimals versus semi-major axis after 100,000 yr of simula- 
tion. The columns of plots from left to right have initial total 
mass in unresolved debris of (sims. la & lb), 0.25 (2a & 
2b), 2.5 (3a & 3b), and 25 M® (4a & 4b). The top row shows 
the control simulations (l-4a) that do not include dynamical 
friction from the unresolved debris. The bottom row shows 
simulations (l-4b) that included dynamical friction. 

The dynamical friction of the background debris signifi- 
cantly damps the eccentricities and causes inward migration 
of protoplanets (the largest planetesimals; see {3.2 \. The 



maximum eccentricity of the protoplanets in simulations 2b 
and 3b is a few times less than the maximum eccentricity 
of the protoplanets in 2a and 3a (Fig. [3|; they are not zero 
due to a few stirred protoplanets at small semi-major axis. 
Most of the protoplanets in 2b and 3b and both protoplan- 
ets in 4b have effectively zero eccentricity. In addition, the 
dynamical friction has an equalizing effect on the planetesi- 
mals, causing them to evolve (grow) as a single population. 
There is little or no identifiable background population of 
small bodies in simulations 1, 2, and 3b compared to a large 
background in sims 1, 2, and 3a. This is especially notice- 
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Figure 4. Semi-major axis versus particle mass after 100,000 yr 
in units of the initial mass of the planetesimals (1.5 X lO^** g). 
Top panels correspond to simulations la, 2a, 3a, and 4a, bottom 
panels correspond to simulations lb, 2b, 3b, and 4b. 

able in simulation 2b, which shows a smooth continuum in 
mass from protoplanets almost 1 x 10^ times the starting 
mass, mo, to planetesimals of mass 2 mo. In the comparison, 
control simulation 2a shows a clear separation between the 
protoplanets (m ~ 1 x 10^ mo) and the background planetes- 
imals (m =1 — 10 mo). The protoplanets in simulation 4a 
have grown large enough that they have agglomerated with 
or scattered all planetesimals within 10-Rh. Simulation 4b 
is quite different: the two large protoplanets have evidently 
migrated inward, sweeping up any particles they encounter 
and leaving no particles with semi-major axis larger than 
0.65 AU. 



3.2 Mixing 

After 10^ yr the protoplanets in the control simulations are 
relatively heterogenous in composition, however, when dy- 
namical friction of the unresolved debris is included, mix- 
ing is suppressed because the dynamical temperature of the 
resolved particles is kept low. Fig. |3] illustrates the compo- 
sitional mixing and migration of the resolved particles af- 
ter 10^ yr. The color coding of the particles in Fig. [s] de- 
picts the composition of the particle with respect to the 
starting composition at its current location. Green particles 
are composed of material predominantly from the particle's 
current location. Blue particles are composed of material 
predominately outside the particle's current location — from 
larger semi-major axis. Red particles are composed of ma- 
terial predominately inside the particle's current location — 
from smaller semi-major axis. Non-primary colors such as 
purple, cyan, or yellow indicate a mixture of the three cate- 
gories listed above. The color of particles that contain mate- 
rial from various locations is determined by mass weighting 



the contributions from each primary color category: green — 
current location; blue — larger semi-major axis than the cur- 
rent location; red — smaller semi-major axis than current lo- 
cation. Using this coloring scheme, a noticable difference in 
the composition of planetesimals between the control sim- 
ulations and those including full feedback from the debris 
becomes apparent by 100,000 yr. 

When dynamical friction from the debris is not in- 
cluded, the protoplanets seem to be evenly mixed. Consider 
the results from simulations la, 2a, and 3a (Fig. [3|. Most 
protoplanets are green, indicating that the majority of their 
composition is from their immediate starting surroundings. 
In each of these simulations there are also a few red and 
a few blue protoplanets, suggesting scattering of both the 
protoplanets and the smaller planetesimals that they have 
accreted has taken place. Simulation 4a also shows compo- 
sitional mixing with three large non-primary colored proto- 
planets (one violet, one purple, and one pink). The colors 
of the protoplanets suggests that they are relatively hetero- 
geneous in composition. In addition, there is no evidence 
of systematic migration of the protoplanets. When dynami- 
cal friction is included, however, there is little indication of 
compositional mixing because all particles have primary col- 
oring. In simulations lb, 2b, 3b, and 4b all particles are green 
or blue. However, there is strong indication of inward migra- 
tion with simulations 2b, 3b, and 4b showing no particles at 
all at large semi-major axis and no red hued particles. The 
blue particles have moved inward from their original loca- 
tion. Most of the mass in these protoplanets originates from 
the protoplanet's original semi-major axis, which is larger 
than the protoplanet's current location. 

Inward migration is an expected outcome from strong 
dynamical friction. There is still significant mass in debris at 
this time (Fig.[7|, thus the protoplanets are all experiencing 
dynamical friction from the background debris. As the pro- 
toplanets migrate inwards to smaller semi-major axis, due 
to circularization from dynamical friction and the accretion 
of debris, the protoplanets gravitationally scatter and stir 
each other because their gravitational spheres of influence 
(10— ISiin) begin to overlap. This keeps eccentricities above 
zero and prevents dynamical friction from shutting off. 

Figures [5] and [6] show composition histograms of the 
most massive protoplanets after 100,000 yr for simulations 
3a and 3b. The composition is binned in semi-major axis. 
The height of a given histogram bar represents the mass frac- 
tion of material from that semi-major axis bin in the pro- 
toplanet. The error bar denotes the current location of the 
protoplanet and the width of the error bar indicates IGRh- 
The histograms are ordered in terms of mass of the proto- 
planet they describe, with the most massive in the upper left 
corner and least massive in the lower right. The protoplan- 
ets from simulation 3a are well mixed, containing, in general, 
roughly equal amounts of material originating both interior 
and exterior to their location (i.e. in most histograms, the 
error bar is in the middle). In the histograms from simula- 
tion 3b, on the other hand, the protoplanets tend to be on 
the left edge, which is consistent with inward migration. 

3.3 Debris Evolution 

In general, the mass of unresolved debris decreases with 
time. The addition of dynamical friction from the unresolved 
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Figure 3. Semi-major axis versus eccentricity of particles after 100,000 yr. The top row shows simulations la, 2a, 3a, and 4a. The bottom 
row shows simulations lb, 2b, 3b, and 4b. The circles with error bars are particles that have reached masses greater than 100 times the 
starting planetesimal mass (mo = 1.5 X lO^** g). The error bars are 10i?H in length, where Ru = (2M/3M,)i/3 a is the mutual Hill radius, 
M is the mass of the protoplanet and M« is the mass of the star. The colors indicate composition with respect to the particle's current 
location. Green particles are made of material predominantly from near the particles' current location, blue particles contain material 
from larger semi-major axis, and red particles contain material from smaller semi-major axis. Non-primary-colored particles represent a 
mixture (mass-weighted) of material from these three categories. 



debris slows the accretion of the debris onto the resolved par- 
ticles. Figure [7]A.-D shows total mass of debris versus time 
for simulations la & b, 2a & b, 3a & b, and 4a & b, respec- 
tively. In general, the mass in background debris drops to 
less than half of the initial value by 1000 yr in simulations 2- 
4a. However, the decline in debris mass is much slower when 
dynamical friction is included. Dynamical friction from the 
background keeps the eccentricities of the planetesimals low, 
and as a result, this slows the accretion of the debris onto 
the planetesimals, meaning that the debris can have a last- 
ing effect on the dynamics of the planetesimals. Recall that 



the mass accreted onto a g iven planetesimal (Eq. 4 in Lein- 



hardt & Richardson 



eivR^27vapj^ , where e is 



2005 1 Sm 

the planetesimal eccentricity, 7? is its radius, a is the semi- 
major axis of its orbit, p is the mass density of the debris 
in the annulus, St is the time since the last debris accretion 



update, and P is the Keplerian period corresponding to a). 
The mass in the debris begins to drop eventually because 
the eccentricities of the protoplanets are small but non-zero 
and as a result dynamical friction from the unresolved debris 
causes inward migration. The protoplanets become crowded 
at small semi-major axis and gravitationally stir each other, 
which in turn increases their eccentricities. The accretion of 
the debris accelerates since the accretion is directly propo- 
tional to the eccentricity. 

The evolution of the debris in simulations la & b differs 
from that seen in the simulations that begin with non-zero 
mass in debris. Surprisingly, even in this case, including dy- 
namical friction of the background debris affects the evolu- 
tion, though in a more subtle way than for simulations that 
begin with debris. In both simulations (la & lb), the mass 
in debris increases (the result of debris-producing collisions) 
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Figure 5. Composition liistograms of the protoplanets — thie most 
massive planetesimals — for simulation 3a after 100,000 yr. The 
center of the red error bar indicates the current location of the 
protoplanet. The extent of the error bars (edge to edge) is lOrjf. 



Figure 7. Total debris mass versus time. The solid lines are the 
control simulations la, 2a, 3a, and 4a (left to right, top to bot- 
tom). The dashed lines are simulations lb, 2b, 3b, and 4b. 




Figure 6. Same as Fig. |5] for simulation 3b, which includes dy- 
namical friction from unresolved debris. In general, the protoplan- 
ets in this simulation are located at the left end of the range of 
their composition histograms, indicating inward migration. 
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Figure 8. Debris image from simulation 3a: mass in unresolved 
debris as a function of time and semi-major axis. The mass in 
debris is shown in logarithmic grey scale from black to white — 
maximum to minimum (zero) mass. 



until a maximum is reached around 1000 yr, at which point 
the mass in debris begins to decrease. The larger planetesi- 
mals are now big enough that fewer and fewer collisions are 
disruptive. The maximum mass in debris reaches just over 
50 and 10 mo for la and lb, respectively. Even though the 
total mass in debris is 0.001 the total mass in planetesimals, 
it is enough to reduce the mass loss from a collision (see ^3.5| 
& Fig.fni). 



In the control simulations the unresolved debris is ef- 
fectively "cleaned-up" at all semi-major axes, leaving effec- 
tively no debris anywhere by the end of the simulations at 
100,000 yr. As an example. Fig. |8] shows the debris image 
of simulation 3a. The mass in unresolved debris is shown in 
grey-scale — black indicates the maximum mass, white indi- 
cates the minimum mass. The sharp transitions in the grey 
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Figure 9. Same as Fig. |8]but for 3b. 

scale are the result of debris producing collisions. In contrast, 
the debris in the dynamical friction simulations is long last- 
ing and is accreted onto protoplanets at late time and only 
at small semi-major axis. Figure [9] shows the debris image 
for simulation 3b. The debris mass drops to zero only for 
the inner most four annuli. The effects of the debris even in 
the inner most annulus last until ~ 4 x 10'* yr. In addition, 
there are effectively no debris-producing collisions. This can 
be seen clearly by looking at the composition histograms of 



the debris annuli (Fig. 10 1. By 100,000 yr, the debris an- 
nuli in simulation 3a (solid lines in Fig. 10 1 are well mixed. 



indicating several debris-producing collisions. However, the 
composition histograms of simulations 3b (dashed lines in 
Fig. 10 1 show no mixing at all. All of the mass in a specified 



annulus is from that annulus. 

Figures [Tl] and [12] shows the time evolution of semi- 
major axis versus eccentricity for particles in simulation 3a 
and 3b, respectively. The eccentricity of the planetesimals in 
simulation 3a is large (0.07) by 10,000 yrs. In comparison, 
dynamical friction keeps the eccentricity of the planetesimals 
low in simulation 3b until there is significant inward migra- 
tion. By 100,000 yrs, there are several protoplanets in 3b 
with significant eccentricity at a < 0.8. This is because in- 
ward migration has caused crowding and gravitational scat- 
tering at small semi-major axis, which has in turn resulted 
in efficient "clean-up" of the debris. 



[Ml] 



Figure 10. Composition histograms of the debris annuli from 
simulation 3a in solid and 3b in dashed after 100,000 yr. The 
semi-major axis of the debris annulus represented in each plot 
from upper left to lower right moves from the inner to the outer 
edge of the simulated annulus. 



mass (M-i dM/dt (X M", where a is a positive number), 
thus, more massive planetesimals grow more quickly than 
less massive planetesimals. In the oligarch phase, the power 
of the relative growth rate, a, changes sign, so more mas- 
sive protoplanets grow more slowly than less massive proto- 
planets, because the velocity dispersion of the planetesimals 
within several Hill spheres of the protoplanet is now depen- 
dent on the protoplanet's mass ( [Kokubo fc Ida|1998| ). 

At early times, (10^ — 10'^ yr), simulations l-4a all have 
Q > 0, indicating onset of runaway growth. After 10* yr, 
the power of the growth rate drops below zero, a ~ —0.6, 
indicating a transition to oligarchic growth. In comparison, 
at early times, (10^ — 10'' yr) simulations 1 and 2b have 
a growth rate that is weakly positively dependent on M, 
a ~ 0.1. At later times {t > 10"* yr), the relative growth 
rate becomes inversely dependent on mass, a ~ —0.5. Simu- 
lations 3b and 4b, which have more massive unresolved de- 
bris initial conditions, do not seem to go through a runaway 
growth regime. The relative growth rate is inversely propor- 
tional to Af even at early time, a ~ —0.4. Both of these 
simulations (3b and 4b) have relatively constant power-law 
slopes for the entire 10^ yr. The growth rate of simulation 
4b increases at very late time due to significant migration 
of the outer protoplanet. 



3.4 Growth Rate 

All of the control simulations go through two growth 
regimes: runaway growth followed by oligarchic growth (see 
solid lines in Fig. |13| |. In contrast, of the simulations that 
contain feedback from the unresolved debris, only the two 
with small amounts of initial debris (lb, 2b) go through two 
phases of growth (black and blue dashed lines, respectively). 
In runaway growth, the relative growth rate increases with 



3.5 Collision Outcome 

The range of outcomes of planetesimal collisions differs dra- 
matically between the control simulations and those that 
include dynamical friction (Fig. |14[ |. In the control cases, 
there is a broad range of collision outcomes from disrup- 
tive debris-producing collisions to perfect merging events 
wherein all the mass from the projectile and target ends 
up in one post-collision remnant. In the control simulations, 
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Figure 11. Snapshots of semi-major axis versus eccentricity for 
simulation 3a. The particle color coding denotes the amount of 
mixing. The size of the particle is proportional to its mass (see 
Fig.|3]caption). The particles with error bars (10 rjj) are at least 
two orders of magnitude more massive than the initial planetesi- 
mal mass. 



100 yr 



100,000 yr 



Figure 12. Same as Fig. |11| for simulation 3b, which includes 
dynamical friction from the unresolved debris. 



6% to 19% of all collisions produce some debris compared 
to to 4% of all collisions in the simulations that include 
dynamical friction. When dynamical friction from the unre- 
solved background is included, any significant background 
mass reduces the eccentricities of the planeteismals (Fig. |3| 
and thus the impact speed of collisions. Therefore, the colli- 
sion outcomes for simulations 2b, 3b, and 4b are effectively 
all perfect merging events. When there is no initial mass in 



Figure 13. Mass of the instantaneous most massive particle as a 
function of time. Solid lines are from the control simulations (la, 
2a, 3a, 4a); dashed lines are results from simulations that contain 
feedback from the unresolved debris (lb, 2b, 3b, 4b). 



debris (simulation lb), the collisions are more violent, but 
the debris that is produced decreases the relative speeds of 
the impactors, reducing the number of debris-producing col- 
lisions compared to the control case (simulation la). 

The trends in collision outcome are consistent with the 
mean impact speeds (Fig. 



15 I. In the control case with no 



initial debris (sim. la) the mean impact speed stays rel- 
atively high during the entire simulation, fluctuating be- 
tween 0.4 and 0.6 iicrit, where Wcrit ~ ^ \J 'ij^ is the speed 
necessary to escape from an object with gravitational bind- 
ing energy equal to the total mass of the colliding system, 
M = Mproj + Mtarg is the Combined mass of the projectile 
and the target, and R = {Rproj + RtargY^^ is the radius of 
the combined mass assuming the target and the projectile 
have the same bulk density. When dynamical friction is in- 



cluded (Fig. 15 sim. lb), the mean impact speed starts at 
about the same value but drops to about 0.2 VctH by lO** 
yrs. As a result of the slower impact speeds, there are fewer 
disruptive collisions in comparison to the control case. The 
trends for the other simulations are similar to those seen in 
simulations la and lb. The control cases (Fig. 15 sim. 2a, 
3a, & 4a) have initial impact speeds around 0.5«crit that drop 
with time, reaching a minimum around 1000 yr, at which 
point the mean impact speed begins to increase again, com- 
ing close to or surpassing the original mean impact speed by 
10*^ yr. The initial mean impact speeds for simulations 2b, 
3b, and 4b are slower and the evolution of the impact speed 
is shallower than their control counter-parts. This explains 
why all of the collisions in these simulations result in perfect 
merging. 

The total number of collisions is similar for all simula- 
tions, ~ A'^ over 10"' yr, where N is the initial number of 
planetesimals in the simulation. If all collisions in a given 
simulation resulted in only one post-collision remnant, the 
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Figure 14. Collision outcome versus time in units of the total 
mass of the system for simulations la, 2a, 3a, and 4a (top row) 
and lb, 2b, 3b, and 4b (bottom row). The outcome equals the 
mass of the largest post-collision remnant divided by the com- 
bined mass of the projectile and target. An outcome of 1 means 
that the projectile and target have merged. An outcome less than 
one means that some debris has been produced by the collision. 
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Figure 15. Histogram of average impact speed for all collision 
events between resolved planetesimals. The histogram bins are 
logarithmic in time. 



number of collisions would equal — 1. Figure [T4] sliows that 
most collisions — even in the control cases — result in perfect 
merging and therefore guaranteed planetesimal growth. 

To highlight the difference in planetesimal evolution in 
both scenarios, Fig. [iGjshows the colhsional growth of one of 
the most massive particles in simulation 3a and 3b. (Fig. 16 
is subtly different from Fig. |13| which shows the instanta- 
neous most-massive particle and does not necessarily depict 
the growth of a single particle.) Every collision involving 
the chosen particle is indicated as a triplet of points: blue 
for the mass of the target, cyan for the mass of the projec- 
tile, and an open black square for the mass of the outcome. 
Effectively all collisions are perfect merging events in both 
simulations, meaning the outcome is a particle whose mass 
equals the mass of the projectile plus the mass of the target. 
In addition, the mass of the target in each collision involving 
the chosen particle is the same as the outcome mass from 
the previous collision, therefore, there is little mass increase 
from accretion of unresolved debris. However, the particle 
in simulation 3a has at least 4 times as many collisions by 
2 X 10'' yr. The majority of the collisions are with projectiles 
that are much smaller than the target. So although the par- 
ticle in simulation 3a has many more collisions, the particle 
in simulation 3b has reached about the same mass by the 
end of the simulation. 



4 DISCUSSION 

The terrestrial planets in our solar system are co-planar and 
on low-eccentricity orbits. In general, it is argued that a 
damping force is required during the last phase of planet 
formation, when protoplanets grow into planets via catas- 
trophic but infrequent collisions with other protoplanets, in 
order to produce a planetary system with low eccentricity 
and inclination. Without a damping mechanism, such as dy- 
namical friction from a significant mass of small particles 
([O'Brien et al. 



2006 Goldreich et al. 20041 or interaction 



with nebular gas ( [Tanaka fc W ard'2004), the infrequent but 
strong gravitational encounters between protoplanets pro- 
duce planets on highly eccentric, excited orbits. 

In this paper we have investigated the hypothesis that 
particle debris damps eccentricities, by conducting a series of 
numerical simulations of the middle phase of planet forma- 
tion (during which planetesimals grow into protoplanets via 
planetesimal-planetesimal collisions) to determine if plan- 
etesimal collisions can produce and/or maintain the neces- 
sary massive background of small debris for the dynamical 
cooling of the large bodies in last phase of planet formation. 
Our simulations suggest that planetesimal collisions do not 
produce enough background material to provide significant 
dynamical friction during the last phase of planet formation. 
In this work we conclude that if the background material is 
responsible for the dynamical cooling of the protoplanets at 
late times, the material had to come from another source be- 
sides planetesimal-planetesimal collisions, or from collisions 
in an even earlier phase of planet formation. 

It is possible that our simplified model for debris ac- 
cretion is missing some important physics: for example our 
model does not include evolution of the unresolved debris, 
nor does it allow for dynamical heating of the background 
debris as a result of dynamical cooling of the large bodies. 
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Figure 16. Mass evolution for one of the largest particles in simulation 3a (left) and 3b (right). Every collision involving the particle is 
depicted in these plots by a triplet of points: one blue point for the target, one cyan point for the projectile, and one open black square 
for the collision outcome. 



In addition, the model does not calculate the gravitational 
focusing of the debris by the protoplanets at late time. We 
model the debris semi-analytically; including heating, colli- 
sional evolution, or migration of the debris may decrease the 
efficiency of the dynamical friction, so the evolution of the 
protoplanets would take more time, but our general con- 
clusions would not likely change. In fact we have already 
run some short simulations where the Vdisp was allowed to 
evolve with the velocity dispersion of the resolved planetesi- 
mals. The results are consistent with the previous statement. 
In the future we may consider increasing the realism of our 
collision model by including evolution of the debris to test 
this conclusion more rigorously. 



5 CONCLUSIONS 

In this paper we have investigated the effect of initial back- 
ground and impact-generated debris on the middle phase 
of terrestrial planet formation (planetesimals evolving into 
protoplanets). We have extended our numerical method 



from Leinhardt & Richardson ( 2005 \ to include composition 
tracking and dynamical friction from unresolved debris. Nu- 
merical simulations using the extended method were com- 
pared to the results of numerical simulations using the orig- 
inal rubble-pile planetesimal collision model, which allows 
for erosion of planetesimals due to collisions but does not 
include dynamical cooling of the unresolved debris. We have 
found that with a mass of initial debris ^ 10% of the mass in 
resolved planetesimals, the planetesimal evolution is quali- 
tatively similar in both cases. If on the other hand there is 
significant initial debris mass 10% of the mass in resolved 
planetesimals) at the beginning of the simulation, the plan- 
etesimal growth modes change when dynamical friction of 
the debris is included. In particular, planetesimals grow con- 
currently when there is no background of smaller resolved 



planetesimals, which is the situation when the dynamical 
friction of the debris is not included. We find no situation 
in which there is enough debris produced from collisions to 
replenish the debris that is accreted onto the planetesimals. 
In addition, the composition of the resulting protoplanets is 
much more homogeneous than the protoplanets in the con- 
trol simulations. When dynamical friction is included, there 
is less initial mixing but significant inward migration at later 
time. 
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